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Many alternative formulations of Einstein's evolution have lately been examined, in an effort 
to discover one which yields slow growth of constraint-violating errors. In this paper, rather 
than directly search for well-behaved formulations, we instead develop analytic tools to discover 
which formulations are particularly ill-behaved. Specifically, we examine the growth of approximate 
(geometric-optics) solutions, studied only in the future domain of dependence of the initial data 
slice (e.g. we study transients). By evaluating the amplification of transients a given formulation 
will produce, we may therefore eliminate from consideration the most pathological formulations 
(e.g. those with numerically-unacceptable amplification). This technique has the potential to pro- 
vide surprisingly tight constraints on the set of formulations one can safely apply. To illustrate the 
application of these techniques to practical examples, we apply our technique to the 2-parameter 
family of evolution equations proposed by Kidder, Scheel, and Teukolsky, focusing in particular on 
flat space (in Rindler coordinates) and Schwarzchild (in Painleve-GuUstrand coordinates). 

PACS numbers: 04.25.Dm, 02.30.Mv, 02.30.Jr 
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I. INTRODUCTION 

Recently developed numerical codes offer the possibil- 
ity of extremely accurate and computationally efficient 
evolutions of Einstein's evolution equations in vacuum 
0. To take full advantage of these new techniques 
to perform an unconstrained evolution of initial data 
and boundary conditions, we must address an unpleas- 
ant fact: many choices for evolution equations and and 
boundary conditions permit ill-behaved, unphysical solu- 
tions (e.g. growing, constraint-violating solutions) near 
physical solutions. 

By way of example, when Kidder, Scheel, and Teukol- 
sky (KST) evolved a single static Schwarzchild hole as a 
test case, they found evidence suggesting that their evolu- 
tion equations and boundary conditions, when linearized 
about a Schwarzchild background, admitted growing, 
constraint- violating eigenmodes 0] 0|- These eigen- 
modes were excited by generic initial data (i.e. round- 
off error); grew to significant magnitude; and were di- 
rectly correlated with the time their code crashed. As 
this example demonstrates, the existence and growth of 
ill-behaved solutions limit the length of time a given nu- 
merical simulation can be trusted - or even run. 

For this reason, some researchers have explored the 
analytic properties of various formulations of Einstein's 
equations [j II [I and boundary conditions 
lEOl |l2l| used in numerical relativity, searching for 
ways to understand and control these undesirable pertur- 
bations. 

In this paper, we discuss one particular type of un- 
desirable perturbation: short-wavelength, transient wave 
packets. [For the purposes of this paper, a transient will 
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be any solution defined in the future domain of depen- 
dence of the initial data slice. Depending on the bound- 
ary conditions, the solution may or may not extend far- 
ther in time, outside the future domain of dependence. 
Inside the future domain of dependence, however, "tran- 
sient solutions" are manifestly independent of boundary 
conditions.] Depending on the evolution equations and 
background spacetime used, these transients can poten- 
tially grow significantly (i.e. by a factor of more than 10^® 
in amplitude). Under these conditions, even roundoff- 
level errors in initial data should produce transients that 
amplify to unit magnitude. Once errors reach unit mag- 
nitude, then guided by the KST results discussed above, 
we expect nonlinear terms in the equations to generically 
cause these errors to grow even more rapidly, followed 
shortly thereafter by complete failure of a numerical sim- 
ulation. In other words, if the formulation and back- 
ground spacetime permit transients to amplify by 10^^, 
we expect numerical simulations of these spacetimes to 
quickly fail. 

In this paper we develop conditions which tell us when 
such dramatic amplification is assured. Specifically, we 
describe how to compute the amplification of certain 
transients for a broad class of partial differential equa- 
tions (first-order symmetric hyperbolic PDEs) that in- 
cludes many formulations of Einstein's equations. If this 
amplification is larger than 10^^, then we know we should 
not evolve this formulation numerically. 



A. Outline of remainder of paper 

In this paper, we analyze the growth of transients. [Re- 
member, in this paper a transient is any solution defined 
in the future domain of dependence of the initial data 
slice.] Rather than study all possible formulations, we 
limit attention to a class of partial differential equations 
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we can analyze in a coherent, systematic fashion: first- 
order symmetric hyperboHc systems. Fm'thermore, be- 
cause we concern ourselves only with stability and the 
growth of small errors, we limit attention to linear per- 
turbations upon some background. Finally, to be able 
to produce concrete predictions, we restrict attention to 
those transients which satisfy the geometric optics ap- 
proximation. 

In Sec. introduce an explicit ray-optics-limit 

solution to first-order symmetric hyperbolic linear sys- 
tems - a class which includes, among its other elements, 
linearizations of certain formulations of Einstein's equa- 
tions. We provide explicit ODEs which determine the 
path (i.e. ray) and amplitude of a geometric-optics so- 
lution, in terms of initial data at the starting point of 
the ray. Then, in Sec. IIIII we introduce wave packets as 
solutions which are confined to a small neighborhood of 
a particular ray. We further define two special classes of 
wave packet - coherent wave packets and prototyptical 
coherent wave packetes - which, because of their simple, 
special structure, are much easier to analyze. Finally, in 
Sec. liVI we introduce and discuss the technique (energy 
norms) we will use to characterize the amplitude of wave 
packets. In particular, we provide an explicit expression 
[Eq. |j2n|)] for the growth rate of energy of a prototypical 
coherent wave packet. 

To demonstrate explicitly how the techniques of the 
previous sections can be applied to produce the growth 
rate of transients, in Sections IVl and IVTI (as well as Ap- 
pendix we describe by way of example how our meth- 
ods can broadly be applied to the two-parameter formu- 
lations that Kidder, Scheel, and Teukolsky (KST) have 
proposed 1]. Specifically, Sections Ivl and IVII will re- 
spectively describe wave packets on flat-space (written in 
Rindler coordinates) and radially propagating transients 
on a Schwarzchild-black-hole background (expressed in 
Painleve-GuUstrand coordinates) . 

Finally, to demonstrate explicitly how expressions for 
the growth rate of transients can be used to filter out par- 
ticularly pathological formulations, in Sections IVIII and 
I Villi we use the results for the growth rates of transients 
obtained in Sections IVland rVll to determine what pairs of 
KST parameters (7 and z) guarantee significant ampli- 
fication of some transient propagating on a Rindler and 
Painleve-GuUstrand background, respectively. 



Guide to the reader 

While the fundamental ideas behind this paper - the 
study of wave packets and the use of their growth rates to 
discover ill-behaved formulations - remains simple, when 
we attempted to perform practical, accurate computa- 
tions, we quickly found the simplicity of this idea masked 
behind large amounts of novel (but necessary) notation. 
We therefore found it difficult to simultaneously satisfy 
the casual reader - who wants only a summary of the 
essential results, and who is still evaluating whether the 



results and the methods used to obtain them are worthy 
of further attention - and the critical reader - who needs 
comprehensive understanding of our methods in order to 
evaluate, duplicate, and (potentially) extend them. We 
have chosen to slant the paper towards the towards the 
critical reader; thus this paper is a comprehensive and 
pedagogical introduction to our techniques. 

While this paper can be consumed in a single read- 
ing, for the reader interested in a brief summary of the 
essential ideas and results, or for anyone making a first 
reading of this paper, the author recommends only read- 
ing the most essential details. First and foremost, the 
reader should understand the scope and significance of 
this paper (i.e. read the abstract and Sec. I). Next, 
the reader should follow the general description of the 
techiques in Sections IIIII and IIVI in detail. Subse- 
quently, the reader should examine our demonstration 
that our techniques indeed give correct results for growth 
rates (cf. the introduction to Section and the sum- 
mary that section's results in Section fV D|l . Finally, to 
understand how these techniques can be used to discrover 
ill-behaved formulations, the reader should examine Sec- 
tions EUl and Eml 

The more critical reader may wish to test and verify 
our computations. This reader should then review Sec- 
tions mil and IIVI again, then work through Sections 
IVl and IVII in detail (returning to the earlier sections for 
reference as necessary) . This reader will also benefit from 
the general approach to KST 2-parameter formulations 
discussed in Appendix IeI 

Finally, the most skeptical readers will want to exam- 
ine the conceptual underpinnings of and justifications for 
our every computation. This reader should simply follow 
the text as presented, but carefully read every footnote 
and appendix as they are mentioned in the text. In par- 
ticular, this reader will want to review our Appendicies 
ini(for a justification of our ray-optics techniques) andlXI 
(for many useful identities used in the previous appendix 
and elsewhere in the paper) as well as Appendix [CI (for 
a more detailed discussion of prototypical coherent wave 
packets, a key element in our computational method). 



B. Connection with prior work 

1. Study a short-time, rather than long time, instability 
mechanism 

First and foremost, we should emphasize that our work 
differs substantially from all previoius work on this sub- 
ject: we very explicitly restrict attention to amplification 
over only a short time (i.e. a light-crossing time). On the 
one hand, unlike other work, because of this restriction, 
our claims - being independent of boundary conditions - 
apply to all boundary conditions. On the other, because 
we forbid ourselves from studying our solutions outside 
the future domain of dependence of the initial data slice 
- even though, in practice, we could draw some elemen- 
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tary conclusions - in this paper we choose not to make 
any claims about how a formulation of Einstein's equa- 
tions will behave at late times (i.e. its late-time stability 
properties). 



2. Study an instability mechanism, not necessarily the 
dominant one 

In other papers which attempt to address the stability 
properties of various formulations of Einstein's equations 
- for example, Lindblom and Scheel (LS) - the authors 
try to (somewhat naturally) an understanding of the 
dominant instability mechanism. Unfortunately, we do 
not fully understand all the dominant instability mech- 
anisms which can occur in generic combinations of evo- 
lution equations, boundary conditions, and background 
spacetimes. Indeed, while some theoretical progress has 
been made towards estimating the dominant instability 
mechanisms (i.e. LS), for generic "reasonable" formula- 
tions (i.e. those which we have not excluded based on 
other known pathologies, such as being weakly hyper- 
bolic), we currently can only reliably determine how ef- 
fective simulations will be by running those simulations. 
And simulations are slow. 

In this paper, instead of studying the dominant insta- 
bility, we study an instability (transients) which we can 
easily understand and rigorously describe. We use this 
instability to discover formulations which are known to 
be troublesome (i.e. which have trouble with transients). 



3. Short-wavelength approximations 

This paper makes extensive use of geometric optics, 
a special class of short-wavelength approximation. Sev- 
eral authors have applied short wavelength techniques to 
study the stability of various formulations of Einstein's 
equations Q, 0, 0- These techniques, however, have 
generally been applied to systems whose coefhcients do 
not vary in space, limiting their validity either to very 
small neighborhoods of generic spacetimes, or to flat 
space. Previous analyzes have thus obtained only a de- 
scription of local plane wave propagation: in other words, 
local dispersion relations. In this paper, with the geomet- 
ric optics approximation, we describe how to glue these 
local solutions together. Such gluing is essential if we are 



In fact, because these solutions are high-frequency solutions, we 
can quite easily determine their interaction with most boundary 
conditions. For example, maximally-dissipative boundary condi- 
tions (i.e. the time derivatives of all ingoing characteristic fields 
are set to zero) imply, in the geometric-optics limit, that all so- 
lutions on ingoing rays will be zero. In particular, that implies 
that, when wave packets reach the boundary, they leave with- 
out reflecting. Other boundary conditions may also be easily 
analyzed. 



to obtain a good approximation to a global solution of the 
PDE and hence a concrete, reliable estimate of the am- 
plification of a transient. In this sense, the present paper 
is the logical extension of work by Shinkai and Yoneda 
(see, e.g., 0), an attempt at converting their analysis to 
precise, specific conditions one can impose which insure 
that transients do not amplify. 

4- Energy norms 

This paper also employs the energy-norm techniques 
introduced by Lindblom and Scheel (LS) 0. Energy 
norms provide a completely generic approach to deter- 
mining the growth rate given a known solution and, 
moreover, can be used to bound the growth of generic 
solutions. While LS choose to apply these techniques 
to study a different class of solution - large-scale solu- 
tions whose growth presently limits their numerical sim- 
ulations - these techniques remain generally applicable. 
We use them to characterize the growth of wave packets. 



II. RAY OPTICS LIMIT OF FIRST-ORDER 
SYMMETRIC HYPERBOLIC SYSTEMS 

In classical electromagnetism, certain short- 
wavelength solutions to Maxwell's equations can be 
approximated by a set of ordinary differential equations 
for independently-propagating rays: a set of equations 
for the path a ray follows, and a set of equations 
which determine how the solution evolves along a given 
ray jl^ . This limit is known as the ray optics (or 
geometric-optics) limit. In this section, we construct 
an analogous limit for arbitrary first-order symmetric 
hyperbolic linear systems. 

A. Definitions 

We study a specific region of 4-dimensional coordi- 
nate space (t,x), on which at each point we have a A^- 
dimensional (real) vector space V of "fields" u £V. 

Inner products: On the space of fields, an inner prod- 
uct is a map from two vectors u, w to a real number with 
certain properties (bilinear, symmetric, and positive- 
definite). The inner product is assumed to be smooth 
relative to the underlying 4-manifold. The canonical in- 
ner product on (i.e. the iV-dimensional dot-product, 
relative to some basis of fields which is defined every- 
where throughout space) is denoted (,), and does not 
vary with space. We can represent any other inner prod- 
uct in terms of the canonical inner product and a map 
S : V ^ V as {u, Sv), where (u, Sv) = {Su, v). 

An operator Q is said to be symmetric relative to the 
inner product generated by S if {u, SQv) = [Qu, Sv) for 
all u, V. In other words, an operator Q is symmetric if it 
is equal to its own conjugate relative to 5, denoted 
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and defined by (u, SQv) = {Q^u, Sv) for all u, v. Equiv- 
alently, the conjugate relative to S may be defined in 
terms of the transpose (i.e. the conjugate relative to 
5"= 1): 

Qt = S-^Q^S . (1) 

First-order symmetric hyperbolic linear systems (FOS- 
HLS): A first-order symmetric hyperbolic linear system 
has the form 

[dt + A''{x, t)da - F{x, t)] uix, t)^0 (2) 

for u(x,t) a smooth function from the underlying 4- 
manifold into the A^-dimensional space of fields, for A" 
and F some (generally space and time dependent^) linear 
operators on that space, and for A" a symmetric operator 
relative to some inner product. 

If more than one inner product makes A symmetric, 
henceforth, when talking about a specific FOSHLS, we 
shall fix one specific (arbitrary) inner product throughout 
the discussion, and therefore some specific S. 

Characteristic fields and speeds: For all 3- vectors Pa, 
A°'pa is symmetric relative to the inner product generated 
by S. It has a set of eigenvalues, eigenspaces, and (for 
each eigenspace) basis eigenvectors, denoted as follows: 

• ujj{t,x,p) are the eigenvalues of A°'pa\ 

• Bj{t,x,p), where j runs from 1 to the number of 
distinct eigenvalues of A"-pa, are the eigenspaces of 
A" Pa, and 

• Vj^a{t, x,p) are some orthonormal basis of eigenvec- 
tors for the space Bj{t,x,p), where a runs from 1 
to the dimension of Bj . 

Because A°'pa is symmetric relative to the inner product 
induced by S, the eigenspaces are orthogonal relative to 
the inner product, and the eigenspaces are complete. Fi- 
nally, at each point (x,p) and for each eigenspace, there 
is a unique projection operator Pj{t,x,p) which satisfies 
PjV = V \i V £ Bj, PjV = if ?; S Bk with k ^ j. 

We require A°-pa and its eigenvalues, eigenspaces, and 
projection operators to vary smoothly over all x"^ and 
Pb in the domain. [We do not demand the eigenvectors 
themselves to be smooth save in the neighborhood of each 
point {x°-,ph): topological constraints may prevent one 
from defining an eigenvector everywhere (i.e. for all pa 
given x°) ■^.] 



^ As a practical matter, we will limit attention in this paper to A" 
and F varying slowly (or not at all) in time; therefore, all time 
dependence in the operators A°- , F, and S may usually be ne- 
glected. For completeness, however, we retain time dependence 
for readers who may wish to apply these techniques to more 
generic systems. 

^ For example, in the first-order representation of the scalar wave 
equation, two of the eigenvectors at each point {x,p) are es- 
sentially vectors transverse to the surface \p\. These cannot be 
extended over the sphere. 



Group velocity and acceleration: We define the group 
velocity V"' and group acceleration Oj a via 

V;{t,X,p} = ■^LOj{t,X,p) , (3) 
aj^a{t,X,p) EE -—UJj{t,X,p). (4) 

We shall make frequent use of an alternative expression 
for the group velocity, Eq. ljA2|l . which is discussed in 
Appendix^ Among other things, Eq. ljA2|l implies 

UJj{x,p) = Vf-pa ■ 

B. Form of ray-optics solution 

We now construct a solution which approximately sat- 
isfies Eq. |5J). Our method works by constructing a set 
of characteristics (i.e. rays), then integrating some am- 
plitude equations along each characteristic (as an ODE) 
to find the amplitudes farther along the ray. 

In this section, we only introduce the results of our 
analysis. In Appendix^ we provide a more comprehen- 
sive justification of our ray-optics approach. 

Ray-optics solution 

Rather than express our solution in terms of the orig- 
inal A^-dimensional variable u, we introduce -|- 1 new 
variables dj ^ and </> and parametrize the original state 
by 

u = ue"^ (5) 

where we further expand u in terms of the eigenvectors 
vi^a of A°'da(t> at each point {t,x): 

u^'^'^di^a{t,x,d(j))vi^a[t,x,d(j)) . (6) 

I a 

[For notational clarity, the arguments t, x, and daip to the 
functions 4>, vi^a, and dia will in the following be usually 
omitted.] 

In terms of these new variables, a ray-optics solution 
is a solution to the following equations, for some fixed j: 

= [dt + V^{x,d(j,)da]<P, (7a) 
= c?;,/3 for I ^ j , and (7b) 

= [dt+ V;da\ d,,a (7C) 
13 ^ ' 

When we substitute solutions to the ray-optics equations 
[Eq. O] back into the original FOSHLS [Eq. ©], as de- 
scribed in detail in Appendix IbI we find the geometric- 
optics solutions are excellent approximate solutions to 
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the original PDE, so long as certain mild conditions con- 
tinue to hold [e.g. the oscillations in (j) remain rapid com- 
pared to any other length or time scale] . 



C. Interpreting the geometric-optics equations 

We introduce the geometric-optics solution precisely 
because it simplifies the PDE ~ in particular, because it 
converts the problem of solving a general PDE [Eq. ^] 
into the problem of solving uncoupled ODEs [Eq. {Tj)]. 
Specifically, these ODEs consist of the the phase equation 
[Eq. (|7a|l ] - which determines the path of the ray leaving 
a point X consistent with initial data for (f> with gradient 
da<j){x) - and the polarization equations [Eqs. Ij7bp and 
((7cl ] - which allow us to propagate the di^a along each 
ray. 

But while these equations are now ODEs, their struc- 
ture is not particularly transparent. In this section, we 
rewrite the phase equation [Eq. (j7a|l ] and the polariza- 
tion equation [Eq. (|7cll ] to better emphasize their prop- 
erties and physical interpretation. 



1 . Path of the ray 

The physical significance of the phase equation [Eq. 
(|7a|l ] becomes much easier to appreciate when it is rewrit- 
ten in first-order form. When we differentiate that ex- 
pression and re-express the result as an equation for 
ka = da4>, we find 



= dtka + V^dbka + [daV^(x,k)]kb 

= dfka + V^dbka - aj^a{x, k) . 



(8) 



[While k does depend on x, because {dk^Vj)kb — the 
last term in the first line does indeed simplify into —cij^a, 
as stated] . Solutions to this PDE may be constructed by 
gluing together solutions to the following pair of coupled 
ODEs for x{t) and k{t): 



dx" 
~dt' 
dka 

dt 



V^ix, k) 

^a.j {x, k) 



(9a) 
(9b) 



By using the definitions of V^" and flaj , we find these are 
precisely Hamilton's equations, using ujj{t,x,k) as the 
Hamiltonian. 

These two equations define the rays (i.e. characteris- 
tics). Given initial data for ka which has fca = ^ai/) in a 
3-dimensional neighborhood of a point, we have a unique 
ray emanating from each point in that neighborhood. So- 
lutions to Eq. (jSJ follow from joining the resulting rays 
emanating from each point in the neighborhood together; 
and solutions for (p [i.e. Eq. H7a|l ] follow by integrating 
the phase out along each ray. 



2. Propagating polarization along ray 

In practice, the polarization equation [Eq. (|7c|l ] is dif- 
ficult to interpret: since it involves spatial derivatives of 
basis vectors, and since we have freedom to choose our 
basis vectors Vj^a arbitrarily within each subspace Bj , we 
cannot transparently disentangle meaningful terms from 
convention-induced effects. 

To constrain the basis and simplify the equation, we 
sometimes choose a basis in the neighborhood of the ray 
of interest which satisfies the no-rotation condition [dis- 
cussed at greater length in Appendix I A 2j : 



{v^,^^,S{dt + A'^da)vj,p;)^Q 



(10) 



where the square brackets denote antisymmetrization 
over a and f3 [i.e. X\^ai3] = (^a/3 — -'^/3a)/2]. The no- 
rotation condition completely constrains the antisym- 
metric part of an operator [i.e. the left side of Eq. (|10|) ]: 
the condition that the basis vectors Vj^a remain orthogo- 
nal constrains that operator's symmetric part; and there- 
fore the basis Vj,a is necessarily completely specified at 
any point along a ray in terms of initial data for the basis. 

Using the no-rotation condition, we find the polariza- 
tion equation becomes the less-arbitrary expression [Ap- 
pendix 



(11) 

) ^3,0 



- \dt + v;da + -dav;{x, d^) ) d,, „ 



Eh 

/5 ^ 



SF+^dtS+^da{SA-) 



where in the above „ is a no-rotation basis. In Sec. 
mil we will use this expression to motivate the definition 
of prototypical coherent wave packets, which have an ex- 
ceedingly simple growth rate. 



D. When do geometric-optics solutions exist? 

Given initial data (say, for ka and dj^a on some initial 
compact region) , we can in practice always find a solution 
to the geometric-optics equations [Eq. {Tj)] valid for some 
small interval St (i.e. by using general PDE existence 
theorems, like Cauchy-Kowaleski). However, for general 
initial data we cannot solve the phase equation [Eq. H7a|l] 
for an arbitrary time T. By way of example, even if we 
find each individual ray [i.e. each solution to Eq. © 
emanating from each initial-data point] emanating from 
our initial data region out to time T, these rays may cross 
before time T, rendering the geometric-optics solution for 
dj_a both singular and inconsistent at the ray-crossing 
point. (A similar problem arises in classical geometric 
optics.) Furthermore, depending on the structure of A"", 
certain rays may not even admit extension to time T (i.e. 
certain rays may be be future-inextendible, precisely like 
rays striking singularities in GR). 
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A proper treatment of these technical compHcations is 
considerably beyond the scope of this paper. In practice, 
we will assume we have chosen initial data so that our 
geometric-optics solution can be evolved to any time T, 
unless it involves transport into a manifest singularity 
(i.e. a singularity of the spacetime used to generate the 
FOSHLS) before time T . Furthermore, we will assume 
the solution is well-behaved - that is, the congruence has 
finite values for ka, V"", Uaj and their first derivatives. 
With a well-behaved solution to the phase equation [Eq. 
(|7a|l ] . we may always find a finite, consistent solution to 
the polarization equation [Eq. H7c|) ] in terms of the initial 
data.^ 



III. DEFINING WAVE PACKETS 

In Sec. ^ we have constructed approximate solu- 
tions to linearized first-order symmetric hyperbolic PDEs 
in the geometric optics limit. These solutions are con- 
structed by integrating ODEs for (and along) rays [Eq. 
{7\. Since each ray evolves independently, we are nat- 
urally led to consider wave packets - that is, ray-optics 
solutions which are nonzero only in a (4-d) neighborhood 
of some (4-space) ray. 

In this section, we outline how wave packets may be 
generally constructed. We also describe the two special 
classes of wave packets, coherent wave packets and pro- 
totypical coherent wave packets (PCWP), which will be 
the focus of discussion henceforth. 



A. Constructing wave packets 

A wave packet that persists for a time T is some so- 
lution to the geometric-optics equations [Eq. {T))] which 
is nonzero only in some small neighborhood of a ray (i.e. 
nonzero only within some coordinate length 5 from the 
central ray). 

From a constructive standpoint, while we can easily 
construct solutions from initial data for ka and dja, we 
have no transparent way, besides solving the equations 
themselves, to determine whether a particular set of ini- 
tial data for ka even generates a congruence which ex- 
ists and remains well-behaved (e.g. 9°Vf, and dakt both 
finite) for time T, let alone whether the specific combi- 
nation of initial data for ka and dja yields a geometric- 
optics solution with support only within a given distance 
i5 from a ray. 

Still, physically we expect we can avoid these technical 
complications. For example, we expect that, for all rays of 



* Since the polarization equation [expressed as Eq. I7cl or as 
Eq. IIII I is linear in the polarization fields dja, it therefore 
admits well-behaved solutions for the evolution of dj^^ along a 
well-behaved ray so long as the linear operators present in that 
equation are well-behaved. 



physical interest, we can extend the central ray of interest 
to time T (i.e. characteristics of physical interest can be 
extended as long as physically necessary) . We expect that 
singular congruences ka can be avoided by proper choice 
of initial ka data (e.g. the ray equations do not require all 
congruences near the ray of interest to diverge or come 
to a focus). And given a well-behaved congruence, we 
expect we can always choose initial data for dja in a 
sufficiently small neighborhood so the solution for dj,a 
is nonzero only within some fixed distance S from the 
central ray. 

Thus, as a proper treatment of these technical compli- 
cations is considerably beyond the scope of this paper, we 
shall henceforth simply assume that a wave packet solu- 
tion can always be constructed about any ray of physical 
interest. 



B. Specialized wave packets I: Coherent wave 
packets 

Since rays propagate independently, one can choose ar- 
bitrary initial data, and in particular arbitrary polariza- 
tion directions w, and still obtain a wave-packet solution. 
Here, w is defined by 

w = u/\u\ \u\ = [{u,Su)]^/^ . (12) 

We prefer to further restrict attention to those wave 
packets which have a single, dominant polarization di- 
rection w present initially (and therefore for all time). In 
other words, we require w vary slowly across the wave 
packet's spatial extent. Wave packets with this property 
we denote coherent wave packets. 



C. Specialized wave packets II: Prototypical 
coherent wave packets (PCWPs) 

While coherent wave packets have a simple polariza- 
tion structure, characterized by some polarization direc- 
tion w, this polarization structure need not necessarily 
have a transparent relationship to the terms present in 
the polarization equation [Eq. (|7c|l : or equivalcntly Eq. 

if we use a no-rotation basis]. Therefore, we de- 
fine prototypical coherent wave packets (PCWPs) as wave 
packets which have at each time their polarization direc- 
tion w equal to one of the eigenvectors /j''^ of the oper- 
ator Oj-. 

O, ^ P, |f + ^S-' [d,S + da{SA^)]j P, (13) 

O,/]^) ^ o,^/j^) (14) 

where /x, running from 1 to the dimension of Bj , indexes 
the eigenvectors of Oj . For simplicity, we assume Oj has 
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a complete set of eigenvectors.^ 

If PCWPs exist, we expect - because of their relation- 
ship to the terms of the polarization equation [Eq. (fTT|l ] 
- the propagation of their polarization will be much eas- 
ier to understand. Most notably, as we will show in the 
next section [Sec. II V| . prototypical coherent wave pack- 
ets have particularly simple expressions for their growth 
rates [i.e. Eq. pHl ]. 

PCWPs will exist as exact solutions to the polarization 
equation [Eq. only in certain special circumstances; 
for example, most of the polarizations to be discussed in 
Sections IVI and IVIl admit exact PCWP solutions. How- 
ever, as demonstrated in more detail in Appendix [CI we 
do not expect the polarization equation to generically 
admit PCWP solutions. 

Nonetheless, as discussed in greater detail in Appendix 
O a PCWP with w = f/^ is a good approximate solu- 
tion to the polarization equation when the eigenvalue Ojv 
of Oj is sufficiently large. Indeed, by rewriting the polar- 
ization equation in the basis fj^\ we can show generic 
coherent wave packets will rapidly converge to a PCWP 
with w = fi'^"'^ for Vo indexing the eigenvalue of Oj with 
largest real part. In other words, based on Eq. (pn|l . 
when coherent wave packets grow quickly, they can al- 
ways be well-described by a PCWP. 



IV. DESCRIBING AND BOUNDING THE 
GROWTH RATE OF WAVE PACKETS 

Since a wave packet is narrow and we care little about 
its precise spatial extent, we commonly characterize the 
wave packet by a single number (e.g. a peak amplitude) 
rather than a generic distribution of polarization over 
space. Unfortunately, the maximum value of the ampli- 
tudes dj,a depend on the spatial extent of the wave packet 
- in other words, it depends on our choice of congruence, 
rather than the central ray itself. 

Because the amplitude function is subject to focusing 
effects (through the term daV""), we choose to describe 
the magnitude of the wave packet by the magnitude of 
its energy norm. Introduced by Lindblom and Scheel 
(LS), the energy norm is an integral quantity analogous 
to energy 0; and, like the energy of a wave packet solu- 
tion to Maxwell's equations, the energy norm will not be 
susceptible to focusing effects. 

In this section, we describe how energy norms can be 
used to characterize the magnitude of wave packets. We 



^ The behavior of the polarization equation when Oj has Jordan 
blocks is straightforward (i.e. we converge to some specific eigen- 
vector in the Jordan block; we obtain no change to the final 
predictions for exponential growth rates; we only add at most a 
polynomial in t to the amplitude functions) but tedious to de- 
scribe in detail. Moreover, in all physically interesting cases we 
have examined, Jordan blocks have not appeared in Oj \ we have 
been able to choose a complete set of basis eigenvectors. 



also obtain special expressions for the growth rates of co- 
herent wave packets [Eq. (|19|l ] generally and prototypical 
coherent wave packets [Eq. (|2()|l ] in particular. 

Also, for completeness, in Appendix IdI we provide an 
explicit, rigorous bound for the growth rate of energy 
which will not be otherwise used in the paper. 



A. Energy norms and the magnitude of 
geometric-optics solutions 

Lindblom and Scheel define the energy norm by way of 
two quadratic functional of a solution u [LS Eqs. (2.3) 
and (2.8)]. When expressed in terms of our notation, 
these functional are terms of our notation are 

e EE {u*,Su) E = Jutfxe. (15) 

Unlike LS, we do not generically have a preferred spatial 
metric; we therefore replace the factor ^ present in LS 
Eq. (2.8) by the more generic fj,.^ 

We may substitute in the expressions appropriate to 
a ray-optics solution to obtain excellent approximations 
to the energy. By way of example, the energy Ej of a 
geometric-optics solution propagating in the jth polar- 
ization may be expressed as 




Sv,.,p) (16) 



a 

where the terms neglected are small in the geometric op- 
tics limit and where the second line holds because by 
construction the basis Vj^a is orthonormal. 



B. Energy norms and the growth rate of wave 
packets 

Following the techniques of Lindblom and Scheel, we 
can use energy norms and conservation-law techniques to 
obtain a general expression for the growth rate of a wave 
packet. 

To follow their program, we must generate a conserva- 
tion law. Define, therefore, an energy current j° [i.e. LS 
Eq. (2.4)] 



Unlike LS, we are not necessarily working with a metric space; 
therefore, we have no preferred measure on the coordinate space 
and therefore allow for an arbitrary, as-yet-undetermined mea- 
sure factor fi. 
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The quantities e and j" obey the conservation-law-form 
equation 

dte + ^-^da inf ) = {u*,SFu) + iFu*,Su) 
+ {u, [dtS + fi-^daifiSA'')] u) 

[i.e. the analogue of LS Eqs. (2.5) and (2.6)]. 

For a wave-packet solution, which is concentrated at 
each time to a small spatial region, the current j" drops 
to zero rapidly, and is in particular zero at the manifold 
boundary. As a result, when we integrate the conserva- 
tion law, we find the energy obeys the equation 

^ = J IJ-d^xiu^SQu) (17a) 

Q = F + S-^[F^S + dtS + ^i-^daifiSA'')]{17h) 

where F"^ is defined so (u, Fv) — {F'^u, v) for all u, v (i.e. 
F^ is the transpose). [In LS, the analogous equations are 
(2.7) and (2.9); in our case, however, we have no surface 
term involving ja because the solution falls off rapidly 
away from the wave packet.] 

We can show Q is symmetric relative to SJ We can 
also show that that Q is closely related to the symmetric 
part of the operator Oj [Eq. (|13|l ]: 



P^QP^ = o, + o] + ^p,v; 



(18) 



C. Energy norms and the growth rate of coherent 
wave packets 

Since coherent wave packets are both localized and 
possess a well-defined polarization direction w, we find 
Eq. (|17(l becomes, for coherent wave packets, 



(19) 



where the right side is evaluated at the location of the 
wave packet at the current instant. 

Because we still need the appropriate polarization di- 
rection w to make use of the above expression - a direc- 
tion we can only obtain from the polarization equation 
[Eq. Hll(l ] - Eq. (|19|l provides only an alternate per- 
spective on the growth of wave packets, not an entirely 
independent approach to the evolution of the amplitude. 



D. Energy norms and the growth rate of PCWPs 

In the special case of a PCWP, however, we do know 
the polarization direction w: it is one of the normalized 



^ Because S and 5^4" are symmetric relative to the canonical inner 
product, so are their derivatives. And if T is symmetric relative 
to the canonical inner product, then S~^T is symmetric relative 
to the inner product generated by S. 



eigenvectors /j^' of the operator Oj [see Sec. IIII C| . In 
this case, we find the energy growth rate for a PCWP 



with 



rf> to be 



1 

Ejfj, dt 



(20) 



[Here, we used Eq. (O in Eq. l(TO|l.] 



V. GEOMETRIC OPTICS LIMIT OF KST: 
RINDLER 

In the previous sections fScctions lTll IIIII and IIV|I . we 

have developed a procedure for computing the evolution 
and amplification of ray-optics solutions in general and 
prototypical coherent wave packet solutions in particu- 
lar. To provide a specific demonstration of these meth- 
ods, we demonstrate how to construct the geometric op- 
tics limit (as described in Section |nj and compute the 
growth rate of wave packets (as described in Sections IIIII 
and II V|) when the first-order hyperbolic system is the 2- 
parameter first-order symmetric hyperbolic system Kid- 
der, Scheel, and Teukolsky introduced (see their Section 
II J) , linearized about a flat-space background in Rindler 
coordinates. 

Our computations in this section proceed as follows. 
First, we review Rindler coordinates and the effects of 
using Rindler coordinates as the background in the lin- 
earized KST equations. We then describe the limited set 
of rays we will study (i.e. rays that propagate only in 
the X direction). Subsequently, we construct the explicit 
form of the polarization equation [Eq. (|7c|) ] for packets 
that propagate only in x. [The analysis simplifies sub- 
stantially because the basis vectors used do not vary with 
x; therefore, the derivatives present in Eq. lf7c|l disap- 
pear.] The analysis of the polarization equation leads us 
directly to an an explicit expression for the growth of en- 
ergy of a coherent wave packet [Eq. H17() ] in general and 
a prototypical coherent wave packet in particular [Eq. 

Finally, to verify our expressions give an accurate de- 
scription of the growth of PCWPs, we compare them 
them against the results of numerical simulations. 



A. Generating the FOSHLS using the bacl4ground 
Rindler space 

Flat space in Rindler coordinates is characterized by 
the metric 



-x^dt^ + {dx"^ -f dy^ + dz^) 



(21) 



for a; > 0. Using this spacetime as a background, we can 
linearize the KST 2-parameter formulation to generate a 
FOSHLS of the form of Eq. |j2Jl - and in particular find 
explicit forms for the operators A"" and F. For example, 
we find that the principal part has the form [KST Eq. 
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(2.59), along with the definition of 5o in KST Eq. (2.10)]: 



dtSgi. 



kij 







(22a) 
(22b) 
(22c) 



As the right-hand sides of these equations are very long, 
we shall not provide them, or an explicit form for F , in 
this paper. 

Using the FOSHLS obtained by linearizing, we can 
proceed generally with any linear analysis, including a 
construction of the geometric-optics limit. 



B. Describing local plane waves by diagonalizing 



The geometric-optics limit is a short- wavelength limit. 
Naturally, then, the first step towards the geometric- 
optics limit is understanding the plane-wave solutions in 
the neighborhood of a point. We find these solutions by 
substituting into Eq. Q the form u (x Uq expi{k-x — Lot); 
assuming k and lo are large, so we may disregard the right 
side; assuming both Uo and A" are locally constant; and 
then solving for Uo and the relationship between ka and 
u. In other words, we find those local-plane-wave solu- 
tions by diagonalizing A'^ka, as discussed in Sec. ^ to 
find eigenvalues Uj and eigenvectors Vj^a, where j indexes 
the resulting eigenvalues and a indexes the degenerate 
eigenvectors for each j . 

Because the principal part is both simple and inde- 
pendent of the two KST parameters (z and 7), we can 
diagonalize it by inspection. For every propagation di- 
rection, the eigenvalues are precisely ijJs{x,k) = s\k\ for 
s = ±1,0. For our purposes, we study only propagation 
in the x direction. Thus, we need only the eigenfields of 
A°-Xa, which are [see KST Eq. (2.61) and also Appendix 

mi 





z.ab 



= 9ab 
= Myab 



U. 



M,ab 

{Pab ± Mxab) 



(23a) 
(23b) 
(23c) 

(23d) 



These expressions may be interpreted as equivalent to 
the basis vectors discussed in Appendix El [see 

Appendix IE 1 cl and in particular Eq. HE6|) ]. 



C. Deriving the polarization and energy equations, 
for propagation in the x direction on the light cone 

In this section, we describe how to construct and an- 
alyze the polarization equation [Eq. (|7c|l ] and energy 
equation [Eq. (|17|) ] for wave packets propagating in the 



X direction. For technical convenience, we limit atten- 
tion to rays which propagate on the light cone ~ in other 
words, which travel on one of the two null curves of the 
metric: 



dx/dt = sx 



for s = ±1. 



1. Essential tool: Diagonalizing P^FPs 

We have the polarization equation [Eq. lf7c|l ] and a ba- 
sis [Eq. (|23|l . or equivalently Eq. (|E6|) ]: the application 
is straightforward. We can, however, substantially sim- 
plify our expression by changing the basis used to expand 
u from Vj^a to the basis of eigenvectors fs^^ of PgFPs, 
defined by the normalized solutions to 

[Equivalently, we may define these eigenvectors in compo- 
nent fashion. For each s, the matrix (fs,a, Fvs^/s) admits 

a complete set of normalized eigenvectors ■ 



/3 



Using these eigenvectors, we regenerate jl 
J2a fs'^aVs,a, which are eigenvalues of PsFPs.] 

These eigenvectors may be classified according to their 
symmetry properties under rotations about the propaga- 
tion axis x: 

• Symmetric-traceless-transverse 2-tensor [basis vec- 
tors correspond to the fields Uy^ and {Uyy — 
U^^)/^/2] One subspace corresponds to the 
2-dimensional space of symmetric-traceless 2- 
dimensional tensors transverse to the propagation 
direction. The operator PsFPs is degenerate in this 
subspace; the single eigenvalue associated with this 
subspace is given by Cs,t, defined by 



(24a) 



Transverse 2-vector [basis vectors correspond to 
the fields U^^ and U^y] Another subspace corre- 
sponds to the 2-dimensional space of 2-dimensional 
vectors transverse to the propagation direction. 
Again, the operator is degenerate on this space. 
The eigenvalue of F in this subspace is given by 
Cs,v for 



Cs,v 



1-^7 
'-l + 2'y 



(24b) 



• 2-scalars [spanned by vectors corresponding to 
the fields U^^ and {U^y + J7|J] Finally, the 2- 
dimensional space of rotational 2-scalars has its de- 
generacy broken by F. For each s, we find two 
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eigenvalues, denoted Cs,si and Cs,s2, with values 

CsM = -s (24c) 
1 + 272 

(s,s2 = ^ , (24d) 



-1 + 27 

These eigenvectors / are linearly independent. Indeed, 
symmetry guarantees that - with the exception of the 
two 2-scalar eigenvectors - most of the eigenvectors are 
mutually orthogonal. 



2. Polarization equation for general geometric-optics 
solutions 



We can apply these eigenvectors to rewrite the polar- 
ization equation [Eq. (|7c|l ] using the basis fs'^K Specifi- 
cally, we define by the expansion ds^a = I]^ ■ 

Noting our basis vectors /i'^' are independent of space 
and time, we find a set of independent equations for the 
Dg^ of the form 



{dt + sxdx)Dsf_, = Csp^sp 



(25) 



This equation, along with the explicit forms for the basis 
vectors /i'^^ , tells us how to evolve arbitrary polarization 
initial data along our congruence. 



3. Energy equation for general geometric-optics solutions 

Similarly, we may rewrite expressions for the energy 
E [Eq. (ifSJ), or Eq. ^] and growth rate E'^dE/dt 

[Eq. lfT7|) ] using the basis For example, we define 

energy of the wave packet by Eq. H16|) . using a measure 
H = — 1 consistent with the flat spatial metric of the 
background. We find, using symmetry properties of the 
eigenvectors to simplify the sum. 



E ^ Jd^'x J2 



(26) 



fie{t.,v} 
j d^x 2Re \d 



The growth rate of energy E ^dE/dt can be obtained in 
two ways: 

1. First, we can explicitly differentiate Eq. H2t)|) . using 
Eq. (|25|l to simplify as necessary. 



2. Alternatively, we can employ the general expression 
for the growth rate of geometric-optics solutions 
[Eq. H19I) ]. [To do so, we express Q in terms of Og 
via Eq. (|18|) . Then we find the following explicit 
expression for Oj by using Eq. IjEQp from Appendix 
IeI which in this case tells us 



PsS-^[dtS + da{SA'')]Ps^sPs 



when we rewrite the results of that expression in 
an operator, rather than component, notation. Fi- 
nally, we employ the basis /i^"*. Because of Eq. 
l(7fjl . we know the eigenvectors /i^'' of PsFPg are 
equivalently eigenvectors of Os-] 

In either case, one concludes 

^= fd'xJ2\D^,n2ReiCs,) + s] (28) 

+ I d'x 2Re [dI,,D,^s2 (/i^^^*, ^/i^')) 
X (C.l + G,.2 + s)] . 

The above equations remain completely generic and ap- 
ply to all ray-optics solutions that propagate along the 
congruence dx/dt = sx. 



4- Energy equation in a special case: PCWPs 

As Eq. H25|) demonstrates, the polarizations do not 
change direction as they propagate. In other words, if 
a wave packet initially has only D^fj, ^ for some spe- 
cific pair of (s, /i), then the wave packet will always have 
Z^s/j 7^ only for that s and /i. Moreover, as noted in the 
discussion surrounding Eq. (|27|) . the basis vectors /i^-* 
used to define the D^^ are eigenvectors of Os ■ Following 
the discussion of Sec. IIII CI we call such a solution a 
prototypical coherent wave packet. 

For a wave packet solution which is confined to the 
(s/i) polarization, we need only one term in each sum to 
find the energy i?sp and growth rate E~i}dEsfi/dt: 



Es^ = / d^x \Ds^\^ 



1 dEs^ 
Esp, dt 



2Re(C^) + s 



(29a) 
(29b) 



[The above expression was obtained directly from Eq. 
(|28|l . Equivalently, we can obtain the same result using 
Eq. (EOl by way of Eq. ^.] 

To be very explicit, we find using Eq. (|24|) the growth 
rates of the tensor (t) and one of the scalar (si) polar- 
izations to be constant, independent of 7 but depending 
on which direction the packet propagates (s = ±1): 



1 dE, 



1 dEs.si 



EsA dt Es^si dt 



(30a) 



We also find the vector (v) and remaining scalar (s2) 
polarizations have a growth rate which varies with 7, 
according to 



— —s 2 1 



Ef,,, dt 



4 + 27 



(27) 



E, 



dt 



-1 + 27 



(30b) 
(30c) 
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D. Comparing growth rate expressions to 
simulations of prototypical coherent wave pulses 

In Eq. H3U|) we tabulated the expected growth rates 
of energy for each possible coherent wave packet. To 
demonstrate that these expressions are indeed correct, 
we compare these predicted growth rates with the results 
of numerical simulations of wave packets propagating on 
a Rindler background. 



1. Specific simulations we ran 

To test the validity of our expressions, we used a ID 
variant of the KST pseudospectral code kindly provided 
by Mark Scheel. He developed this code to study the 
linearized KST equations on a Rindler background (e.g. 
to produce the results shown in Lindblom and Scheel Sec. 
IV A ,2|). 

We ran this code at a fixed, high resolution (512 collo- 
cation points in the x direction) on a computational do- 
main X e [0.01, 1] with various wave-packet initial data. 
Specifically, we used a wave packet profile proportional 
to 

W{x) = Acos{2ttx/\) exp [-{x - Xc)^/cr^] (31) 

with A = 10"^ Xc = 0.55, a = 0.1, and A 0.01. The 
precise initial data used depended on the polarization wc 
wanted: 

• Tensor When we wanted a tensor polarization, we 
used initial data for a single left-propagating 2- 
tensor component: U~y = W, with all other charac- 
teristic fields zero. In other words, we used initial 
data Pyz = Mxyz — W{x)/2 with all other fields 
zero. 

• Vector When we wanted a vector polarization, we 
used initial data for a single left-propagating 2- 
vector component: = W, with all other char- 
acteristic fields zero. In other words, we used initial 
data Pxz = Mxxz = W/2 with all other components 
zero. 

• Scalar 1 (si) When we wanted to excite the left- 
propagating si polarization, we used initial data 
P.x = Mxxx = W/2. 

• Scalar 2 (s2) After some algebra, one can demon- 
strate that to excite the s2 polarization, we should 
use initial data Pyy — M^yy — W/A and Pyy = 

Mxyy = -W/A. 

To avoid the influences of boundaries, we only studied 
the results of the simulations out to a time t ~ 0.1. 




-4 -2 2 4 

Y 



FIG. 1: The two solid curves show the theoretically-predicted 
growth rates for the "vector" [v, Eq. I|3flb^ ] and one of the 
scalar [s2, Eq. ISOcU l polarizations, when those polarizations 
propagate to the left (s = — 1). The circles show the results for 
numerical simulations of the vector wave packet; the crosses 
show the results for wave packets in the s2 polarization. Both 
predictions agree very well with simulations. 



2. Results 

For each polarization (i, v, si, and s2), we found that 
wave packets remained in the initial polarization, with 
little contamination from other fields. For example, when 
exciting the tensor polarization, we found all fields other 
than Uxy remained small. 

The wave packets' energy grew exponentially, with 
growth rates that agreed excellently with Eq. (|30|l . For 
example, the polarizations si and t both had growth rates 
consistent with unity to a part in a thousand. Our ex- 
pressions for the growth rates for s2 and v also agreed 
well with the results of numerical simulations, as shown 
in Fig. ^for left-propagating pulses (s — —1). 



VI. GEOMETRIC OPTICS LIMIT OF KST: PG 

In this section, we study another example of the ge- 
ometric optics formalism: the propagation of radially- 
propagating wave packets evolving according to the 
KST 2-parameter formulation of evolution equations, lin- 
earized about a Painleve-GuUstrand background. 

Our analysis follows the same course as the Rindler 
case addressed in Sec. We first review Painleve- 
GuUstrand coordinates and the effects of using these co- 
ordinates as the background in the linearized KST equa- 
tions. Subsequently, we construct the explicit form of the 
polarization and energy equations [Eqs. ifTTl) and l(T7|l ] 
for packets that propagate radially on the light cone. Fi- 
nally, in a departure from the Rindler pattern, we also 
add an analysis of the "zero-speed" modes that propagate 
against the shift vector. 
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A. Generating the FOSHLS using a background 
Painleve-Gullstrand space 

A Schwarzchild hole in Painleve-Gullstrand coordi- 
nates is characterized by the metric 



dr 



-dt I + r^dn^ . (32) 



We shall use this metric in cartesian spatial coordinates 
[i.e. z — rcosO, x = rsin0cos0, y — r sin sin;/)] as the 
background spacetime in the KST equations. Linearizing 
about this background, we obtain the explicit FOSHLS 
we study in the remainder of this section. 

As before, we shall not provide the very complicated 
derivative-free terms (i.e. i^) explicitly in this paper. The 
principal part, however, remains simple by design; in this 
case, we have [KST Eq. (2.59), along with the definition 
of So in KST Eq. (2.10)]: 



{dt - /3'^9a)5y 
{dt - fi''da)Mk^ + ^kP^J 



(33a) 
(33b) 
(33c) 



with = 1/27^ 



B. Local plane waves and diagonalizing A°'fa 

As discussed generally in Sec. and by way of a 
Rindler example in Sec. IV Bl to understand how wave 
packets propagate radially we must first understand how 
local plane waves propagate radially, which in turn re- 
quires we diagonalize A'^Va- The basis vectors and eigen- 
values are addressed in detail and in a more general 
setting in Appendix IE 1 cl In brief, the eigenvalues are 
ijJs{^, k) = s\k\ — P°'ka with s — ±1, and the eigenvec- 
tors correspond directly to the Rindler results [i.e. Eq. 
(123), with X —f r; the similarity exists because we can 
use symmetry without loss of generality to demand the 
ray propagate radially in the x direction, along r — x]. 



C. Deriving the polarization and energy equations, 
for radial propagation on the light cone 

Almost half (12 of the 30 characteristic fields) naturally 
are associated with wave packets that propagate at the 
speed of light of the background spacetime (i.e. s = ±1). 
In other words, they propagate on characteristics that 
correspond to null curves of the Painleve-Gullstrand met- 
ric [Eq. (|32|l ]. For radially propagating characteristics, 
that means 



yr 




(34) 

(35) 




FIG. 2: Examples of the three types of radially propagat- 
ing rays of the KST 2-parameter system linearized about a 
Painleve-Gullstrand background. The solid lines show rays 
propagating inward at the speed of light The dotted 

lines show rays propagating "outward" at the speed of light 
(V^). Finally, the dashed curve shows the rays that propa- 
gate inside the light cone (at speed Vr)- The quantities Vr 
are defined in Eq. (1351 . 



with s = ±1. The resulting null curve structure is shown 
in Fig. El 

Because both this case and the Rindler case discussed 
in Sec. IV Gl possess rotational symmetry about the prop- 
agation axis, the equations governing these two cases 
prove exceedingly similar. The analysis follows the same 
course. 



1. Essential tool: Diagonalizing PsFPa with s = ±1 

As in the Rindler case, we will rewrite the polariza- 
tion and energy equations by using eigenvectors /i'^'' of 
PsFPs- Because we again have rotational symmetry 
about the propagation direction, we can again decom- 
pose the eigenvectors into a set of two scalars (si and 
s2), a 2-vector v, and a symmetric-traceless-2-tensor t. 
The eigenvalues may be expressed using 



c 



S,fJ, 



2r 



,3/2 



where the Csp are defined by 



Cs.sl 
Cs,s2 



7 (33 -I- 91z H- 24z^) 

2 + ^ ~ 4(l + 3z)(l-27) 
3 — 3z — 57 



1-27 



Cs.t = 1 



(36) 

(37a) 
(37b) 

(37c) 
(37d) 



The eigenspaces are, by symmetry, spanned by precisely 
the same fields as in the Rindler case. In particular, as 
in the Rindler case the eigenvectors do not change as we 
move along a ray. 
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2. Polarization equation for s — ±1 

For polarizations which propagate radiaUy on the hght 
cone (i.e. s = ±1), the polarization equation [Eq. (f7c|l ] 
can be written as 



= [dt+V,^dr]ds.a+S^ 

r 



(38) 



where we make use of Eqs. (jE7|l and IjESp to simplify the 
right side, and where we observe dar"" = 2/r. 

As in the Rindler case, we may expand the ampli- 
tude u = '^fj^Ds^ifs^^ in terms of the basis fi'^\ and 
thereby arrive at polarization propagation equation pre- 
cisely analogous to the Rindler result [compare with Eq. 



dt + 




(39) 



These equations may be integrated to describe the evo- 
lution of polarization along any individual radial ray. 



3. Energy equation for s = ±1 

Because symmetry guarantees a close similarity be- 
tween this Painleve-GuUstrand case and the Rindler case, 
we find the energy E oi a. geometric-optics-limit solution 
propagating on the light cone radially inward (s = —1) 
or outward (s = -|-1) can be expressed with precisely the 
same expression we used in the Rindler case: Eq. (|26ll . 
[In this case, we again use a measure ^ — 1 compati- 
ble with the background flat spatial cartesian-coordinate 
metric] 

The rate of change of this energy, dE/dt, can be ob- 
tained in two ways. On the one hand, we can directly 
form i?, convert to spherical coordinates, differentiate 
the resulting expression for dE/dt, and use Eq. 
On the other hand, we can find dE/dt using the general 
expression of Eq. (|17|l , an expression we simplify by us- 
ing i) the relation between Q and Oj given in Eq. (|18|l . 

ii) the basis /i'^'' of eigenvectors of PsFPg, and iii) the 
expression [obtained from Eq. (|E9|) and converted from 
a component to an operator expression] 

[dtS + dais A-)] - -^/i^j^P^ ■ (40) 

In either case, we conclude 



,2Re(Gp)-3 



^ = / d^xy\Ds,\ ^- 

+ J d'x 2Re [dI,,D,^s2 (/i^^^*, SI 

Csl + Cs,s2 - 3 



(41) 



In particular, for prototypical coherent wave packets - 
that is, wave packets where s and \i are the same every- 
where in the packet - we can express the growth rate of 
the energy Eg^, of the wave packet as 



1 dEs, 
Eg^ dt 



2Re(C. 



V2r3/2 

where r is the current location of the packet. 



(42) 



D. Deriving the polarization and energy equations, 
for radial propagation against the shift vector 

The remaining 18 fields propagate inward against the 
shift vector, at speed Vo — —y^2/r. 

We shall not follow the same pattern we used to ad- 
dress propagation on the light cone [on a Rindler back- 
ground in Sec. IV CI and on a Painleve-GuUstrand back- 
groimd in Sec. IVI Cj . In those sections, we provided 
extensive discussion and background - the explicit form 
of the polarization equation; a modified form of the polar- 
ization equation in an alternative basis; explicit expres- 
sions for the growth rate of energy general geometric- 
optics solutions; explicit demonstration that PCWP so- 
lutions existed - before finally recovering the growth rate 
of PCWPs. Instead, for pedagogical and other reasons 
[see Sec. IVI D 3| . we shall take a briefer, more practical 
approach better suited to extracting precisely the infor- 
mation needed to decide when some coherent wave packet 
can amplify a significant amount within the future do- 
main of dependence. 

Specifically, following the arguments at the end of Sec. 
nil CI we expect that - whether or not PCWPs exist as 
exact solutions to the polarization equation - when the 
largest eigenvalue Oou of Oo is particularly large, a generic 
coherent wave packet will rapidly converge to a PCWP 
with w = fi^\ In other words, we expect that when 
the growth rates are large, the growth rate of generic 
coherent wave packets can be obtained by finding the 
largest value of dE/dt/ E for PCWPs [i.e. the maximum 
of Eq. ipn|) over ^] . 

In short, we continue to evaluate Eq. H2U|) to get 
growth rates, though now we trust the results only when 
the growth rates are large. 



1. Growth rate of PCWPs 

To evaluate the growth rate of PCWPs, we must diag- 
onalize Oo - 

Oo = Po |P + \s-^ [dtS + daiSA'^)]^ Po 

However, from Eq. (|40|l we know the term in square 
brackets is diagonal. Therefore, diagonalizing Oo to ob- 
tain eigenvalues Oo^i and eigenvectors /j'^'' is equivalent to 
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diagonalizing PqFPo for eigenvalues C,o^J. and eigenvectors 
/i''-* . The eigenvalues of the two operators are related 

by 



Helicity-2 a 4-dimensional space of symmetric- 
traceless-2-tensors ("helicity-2" states), with 
doubly-degenerate eigenvalues given by 



Con 



2\/2r3/2 



(43) 



We shall express the eigenvalues Co^ of PqFPo in terms 
of the dimensionless rescaled quantities , defined im- 
plicitly by 



Com - X V2/r3/2 



(44) 



Substituting Eq. (|43|l into the general expression for the 
growth rate of PCWPs [Eq. wc find that a PCWP 

in the polarization /i will have energy grow at rate 



1 dEo, 

Eon dt 



2Re(L^ 



V2_ 

j,3/2 



(45) 



where r is the instantaneous location of the packet. 



2. Essential tool: Diagonalizing PoFPo 

To obtain explicit growth rate expressions using Eq. 
(|45|l . we need the eigenvalues of PqFPo, expressed ac- 
cording to Eq. 

As in the previous two cases, the eigenspaces of PqFPo 
may be decomposed into distinct classes, depending on 
their symmetry properties of rotation about the propa- 
gation axis. These spaces are as follows: 

• Helicity-0 a 4-dimensional space of rotational 
scalars ( "helicity-0" states) , with eigenvalues given 
by Eq. with 



-l + 3z-|-18i2±y3yiT 



J(-30 



4(H-3z) 
f 19r/+12z 



2z) 



(46a) 
(46b) 
(46c) 



Here, we use 77 = —2/(7 — 1/2) and Yi given ac- 
cording to an expression listed below [Eq. (|47|l ]. 

• Helicity-1 an 8-dimensional space of rotational 
2- vectors ("helicity-1" states), with doubly- 
degenerate eigenvalues given by 



Lyi — —2 



Lv2 
Lv3.v4 



= 5(1 + 6-*) 

8 32(1 + Sz) 



± 



32(1 + 3i) 



(46d) 
(46e) 



(46f) 



Again, we use rj = —2/(7 — 1/2). The expression 
for Y2 is given below [Eq. 



Lti = -(1 + 22) 
Lt2 = 2 + 3z 



(46g) 
(46h) 



• Helicity-3 and finally a 2-dimensional space of 
helicity-3 states, with eigenvalue 

L3 - 3(1 + z) (46i) 

In the above discussion, Yi^2 are defined by 

Yi = (l + 3z)(-5 + 5z + 24z2 + 36z3) (47) 
Y2 = 1296(1 + 3z)2 + 7^2(13 + 83z + 84z2)2 

-247^(1 + 3z)(89 + 199z + 132^2) (48) 

and we use the shorthand 77 = ^2/(7 — 1/2). 

3. Aside: Why can't we follow the previous pattern? 

Unlike all cases previously discussed, a handful of the 
eigenvectors depend weakly on position. As a result, the 
use of a basis which diagonalizes Oo does not offer as 
dramatic a simplification as it did in our earlier analyzes 
of the polarization equation [Sections IV CI and IVI C| . To 
be explicit, if we rewrite the polarization equation in the 
basis fo'^^ in the fashion of those earlier analyzes, we 
obtain [see Eq. (|Cl|l ] 

^ Do,M^, - (dt + V^da + \daV: - o,^ (49) 

for Af^j, some nonzero, position-dependent matrix cou- 
pling the various I?o/i- 



VII. TRANSIENTS AND LIMITATIONS ON 
NUMERICAL SIMULATIONS: RINDLER 

In earlier sections, we developed ~ in general [Sec- 
tions ITU UTTl and IIV| and for specific examples [e.g.. Sec. 
IVl analyzes propagation of transients according to the 
KST 2-parameter formulation of Einstein's equations, 
linearized about a Rindler background] - tools to an- 
alyze the growth of special (i.e. prototypical coherent 
wave packet) geometric-optics-limit transient solutions. 
In this section, we demonstrate how these tools can be 
used to discover when a particular formulation of Ein- 
stein's equations [here, some specific member of the KST 
2-parameter system] which is linearized about a specific 
background [here, flat space in Rindler coordinates] ad- 
mits some massively-amplified transient solution. 

Specifically, in this section we apply the general tools 
developed in an earlier section [Sec. to determine the 
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B. Amplification expected 

For each ray that propagates on the light cone 
(dx/dt = ±a;) within the future domain of dependence, 
and for each polarization on that ray, we can compute 
the amphfication in energy. If Rg^^ = Ej^dEsf^/dt [see 
Eq. 13U|l ] , we can express the ratio of energy of the wave 
packet when it exits the future domain of dependence at 
time iout to the initial energy at time t — Q as 

As^J.{xo) = Esf^itout)/Esf^{0) = exp (iout^s^) ■ 



FIG. 3: The shaded region is the future domain of depen- 
dence of the region x G [0.01, 1] for the KST 2-parameter 
formulation of evolution equations linearized about a Rindler 
background. Transients are any solutions which are defined 
in this region. We study all the prototypical coherent wave 
packets which propagate on the light cone (i.e. according to 
dx/dt — ±x). 



largest possible amplification of a prototypical coherent 
wave packet while it remains within the future domain 
of dependence of some initial data slice. In Sec. IVll Al 
we describe the initial data slice we chose and the sub- 
set of transient solutions we studied. In Sec. IVII Bl we 
apply the tools developed in an earher section [Sec. 
to determine the amplification of each transient. We also 
find an expression for the largest possible amount a tran- 
sient can amplify. Finally, in Sec. IVII Bl we invert our 
expression to determine which pairs of KST parameters 
(z, 7) admit transients that amplify in energy by more 
than 10'^^ (i.e. in amplitude by more than 10^^). 



We have explicit expressions for Rs^,', we can compute 
touti^o, s) for each initial point Xq and for each propaga- 
tion orientation (i.e. for each s); and we therefore can 
maximize Asfj,{xo) over all possible choices of initial lo- 
cation (xo), propagation direction (s), and polarization 
(/i) to find the largest possible ratio A of initial to final 
prototypical coherent wave packet energy. 

In fact, because for each polarization of prototypical 
coherent wave packet, the growth rates of energy is in- 
dependent of time and space, the largest amplifications 
possible always occur along the longest-lived rays - in 
other words, along the two bounding rays x+ and a;_, 
which both extend to tout = ?max- Therefore, we con- 
clude that, while within the future domain of dependence 
of the slice x e [0.01, 1], the largest amount the energy 
of any prototypical coherent wave packet can amplify is 
given by the factor 

A = exp(rmax^Rind) (51) 

where i?Rind is given by 



A. Transients studied 

We limit attention to the future domain of dependence 
of the initial-data slice x G [0.01,1] at i = 0. Since 
the KST 2-parameter formulation has fields which prop- 
agate at (but no faster than) the speed of light, the 
future domain of dependence of this slice is precisely 
what we would obtain using Einstein's equations: a re- 
gion bounded by the two curves X- = 0.01 exp i and 
a;+ = exp(— i). This region is shown in Fig. |21 The 
future domain of dependence extends to time 

T„,ax = lnlO, (50) 

at which point the two bounding curves curves intersect. 

Geometric-optics solutions are defined on rays [i.e. so- 
lutions to Eq. |(SJ)]. While three classes of rays ex- 
ist in this region - those ingoing at the speed of light 
[dx/dt = — x); those outgoing at the speed of hght 
{dx/dt ~ +x); and those which have fixed coordinate 
position - we for simplicity chose to study only the am- 
plification of transients that propagate on the light cone. 



R 



Rind 



maxi? 
max ( 1, 



max[2Re(Cs„) -I- s] 



(52) 



, 1 + 7 

'-I + 27 



. 1 + 27^ 

'-I + 27 



KST formulations which definitely possess some ill-behaved 
transient solution when linearized about Rindler 



Finally, we can invert Eq. H51|l to find those combi- 
nations of KST parameters (z, 7) which permit some 
transient (in particular, some prototypical coherent wave 
packet) to increase in energy by more than a factor 10'^^ 
(i.e. 10^^ in amplitude). The condition may be expressed 
as either A > 10"^^ or, equivalently, as i?Rind > 32. The 
function i?Rind is shown in Fig. 01 along with the line 

-RRind — 32. 

Therefore, we know that some transient can amplify 
in energy by more than 10^^ if i) 7 > (33 + V949)/4, 
n) 7 < -(31 + \/T077)/4, or iii) 7 > 29/64 and 7 < 
(33- \/949)/4. 
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FIG. 4: The solid curve is the theoretical prediction for the 
largest growth rate of wave packets that propagate on the 
light cone [Eq. 1521 1. The horizontal line is the value 32. 
According to arguments made in Sec. IVll Bl those 7 which 
have -RRind > 32 have some prototypical coherent wave packet 
which, in the future domain of dependence, amplifies in energy 
by more than 10'^^. 



C. Relevance of our computation to numerical 
simulations 

We have demonstrated that the KST 2-parameter for- 
mulation of Einstein's equations always admits, at any in- 
stant, prototypical coherent wave packet solutions which 
grow exponentiaUy in time. Generically, we expect that 
at each instant (including in the initial data) these so- 
lutions are excited by errors in the numerical simulation 
(e.g. truncation and roundoff). They then propagate and 
grow; eventually, they reach the computational bound- 
ary. 

Our calculations above describes the largest amount 
any such wave packet solution could possibly grow by 
the time it reaches the computational boundary. If that 
amplification factor is sufficiently large that the wave 
packets reach "unit" amplitude (i.e. whatever magni- 
tude is needed to couple to nonlinear terms strongly), 
here conservatively assumed to be 10^^, then we expect 
any simulation using that particular combination of KST 
parameters will quickly crash. 



Aside: What happens to PCWPs at late times? 

Eventually, the wave packets excited by numerical er- 
rors will reach the computational boundary. What hap- 
pens afterward depends strongly on the precise details of 
the boundary conditions. 

For example, maximally-dissipative boundary condi- 
tions (i.e. the time derivatives of all outgoing char- 
acteristic fields are set to zero) will allow the wave 
packet to leave the computational domain entirely (with 
some small amount of reflection that goes to zero in the 
geometric-optics limit). In this case, at late times no 



transient will ever amplify by more than the amount de- 
scribed above (in Sec. IVll Ij() . 

On the other hand, other choices for boundary condi- 
tions could cause wave packets to reflect back in to the 
computational domain. In these circumstances, the out- 
comes are far more varied ~ at late times, the wave packet 
could potentially grow, could decay to zero, or could en- 
ter a repetitive cycle where on average its amplitude is 
constant.* 

Therefore, without some more specific proposal for 
boundary conditions, we cannot make useful statements 
regarding the late-time development of this instability 
process - or, in other words, we cannot study the growth 
of coherent wave packets for more than a light crossing 
time. 



VIII. TRANSIENTS AND LIMITATIONS ON 
NUMERICAL SIMULATIONS: PG 

In this section, we provide another example of how 
tools developed earlier for the analysis of prototypical co- 
herent wave packets - in general [Sections ITU IIIII and lIV| 
and for specific examples [e.g.. Sec. IVll anal vzes propaga- 
tion of transients according to the the KST 2-parameter 
formulation of Einstein's equations, linearized about a 
Painleve-Gullstrand background] - can be applied to dis- 
cover which formulations of Einstein's equations permit 
ill-behaved transients. 

Specifically, in this section we study the propaga- 
tion coherent wave packets in the 2-parameter KST 
form of Einstein's evolution equations, linearized about 
Schwarzchild written in Painleve-Gullstrand (PG) coor- 
dinates. The theory needed to understand the propa- 
gation and growth of radially-propagating coherent wave 
packets has been developed in an earlier section [Sec. IVI] . 
We apply our techniques to a handful of coherent wave 
packet transient solutions, to discover conditions on the 
two KST parameters (z, 7) which permit amplification of 
those transients' energy by a factor l/e^ for = 10""'^^. 

To provide concrete examples of estimates, we assume 
the initial data slice contains the region r S [2, 10]. So 
any influence from boundary conditions cannot muddle 
our computations, we limit attention to coherent which 
are defined in the future domain of dependence of that 
slice. 



In fact, in this particular case, we expect that if a wave packet 
with growth rate l/r reflects, then symmetry and the structure 
of the Rindler growth rates [i.e. Eq. lijOI I insures that the re- 
flected ray has growth rate — 1/t. Therefore, on average, the 
wave packet has a zero growth rate. 
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FIG. 5: The shaded region is the future domain of dependence 
of the region r € [2, 10] for the KST 2-parameter formulation 
of evolution equations. Transients are any solutions which are 
defined in this region. For reasons emphasized in the text, the 
rays that last for the longest coordinate time prove particu- 
larly helpful. These rays are the left and right boundaries (i.e. 
the horizon and a ray that propagates inward at the speed of 
light from r = 10) and one ray propagating against the shift 
vector which emanates from their intersection. 



A. Transients studied 

We limit attention to the future domain of dependence 
of the region r € [2, 10] at t = 0. Since the KST 2- 
parameter formulation has fields which propagate at (but 
no faster than) the speed of light, the future domain of 
dependence of this slice is precisely what we would obtain 
using Einstein's equations: the region shown in Fig. El In 
particular, the future domain of dependence is bounded 
on the left by the generators of the horizon (trapped at 
r = 2) and on the right by rays travelling inward at the 
speed of light. This ingoing ray reaches r — 2 aX the 
endpoint of the future domain of dependence, at time 

t : 



imax defined by 



10 



dr 



12 1 + ^J2/r 
= 4[3-\/5 + csch"\2)] 



(53) 



4.98 



In our future domain of dependence, we have three 
classes of solutions to the ray-propagation equation [Eq. 
©]: those ingoing at the speed of Hght {V- = — 1 — 
■y/2/r); those ingoing with the shift {Vo = —yj2/r); and 



congruence and in the polarization /i, the energy at the 
time tout{foi s) the ray exits the future domain of depen- 
dence is larger than the initial energy by a factor 



As^^iro) = Esf,{tout)/E,^{0) (54a) 

/■tout rT-a^t 

InAp(ro) - / dtRs^= —Rs^ (54b) 



those outgoing (V+ 



1 



^) [Eq. jSHJ] 



We then search over all Tq, over all propagation direc- 
tions s, and over all polarizations fi to find the largest 
amplification factor A. 

In fact, as in the Rindler case, we immediately know 
which rays produce the largest possible amplification, so 
we can perform the maximization by inspection. 

• Outgoing at light speed: Since the amplification 
of energy increases as r gets smaller {dE/dt/E cx 
l/r"^/^) and with the duration of the ray in time, 
manifestly the generator of the horizon - with both 
the longest duration and the smallest r of all out- 
going rays -will provide the largest possible ampli- 
fication. 

Since the ray of interest has fixed radial location 
r = 2, we find C+,fi is constant for all polarizations. 
Thus, the energy of a prototypical coherent wave 
packet in polarization fi increases by a factor A+p, 
for = exp(t,„axC-i'A')- other words, 

h-iAf, ^ [2Re(C+,p)-3] [3-\/5 + csch~^(2)] (55) 
« 1.245 [2Re(C+,^) - 3] 

[The values for each are given in Eq. (|37|l .] 

• Ingoing at light speed: The longest ray - namely, 
the right boundary of the future domain of de- 
pendence - permits the greatest possible amplifi- 
cations. Thus, among all possible ingoing rays, the 
largest amplification factor for the polarization /i is 
given by A-^: 

r . 1 ln5 - 2csch~^(2) , , 
ln^_^ = [2Re(C-,^) - 3] ^ (56) 

« 0.323 [2Re(C-,^) - 3] 

[The values for each C-fj. are given in Eq. 137|l .] 

Note that since C-^i — C+fj.^ the outgoing transients 
trapped on the horizon grow more than the ingoing 
ones over the same time interval^. 



B. Amplification conditions 



• Ingoing with lapse: The amplification of energy in- 
creases both with ray length and with proximity 



For each of the three classes of rays (s = ±1, 0) propa- 
gating radially in the future domain of dependence [Fig. 
|5] and for each polarization on that ray, we can compute 
the amplification in energy using Rg^ = Ej^dEg^/ dt [see 
Eqs. H42|l and H45|l ]. Specifically, for a wave packet start- 
ing at r = To at time t = 0, propagating in the s-type 



This should be expected: the ingoing and outgoing wave packets 
have similar growth rates at any given radius; we limit attention 
to rays which persist for a fixed time; and the outgoing modes 
we study remain closer to the horizon, where the growth rate is 
larger. 
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to r = (since growth rates go as 1/r^/^). Thus, 
the longest ray propagating at this speed contained 
in the future domain of dependence gives the best 
chances. That ray starts with r ~ r^, with 
defined so the ray terminates at the horizon at 

i — imax • 



40 



20 



(4 + 3W)' 



n 1/3 



(57) 



Thus, we find the largest possible amplification 
among those polarizations that have s = to be 
given by A^)^,: 



2Re(L^) - - 



2Re(i^) - 



dr \/2 
ln(rL/2) (58) 



[The values for each are given in Eq. (|46|l .] 



C. Results: Some KST parameters which have 
transients which amplify by 10''^ 
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FIG. 6: The shaded region indicates those KST param- 
eters which produce some radially-propagating prototypical 
coherent wave packets which amplify their energy by greater 
than 10"^^ within the future domain of dependence of the slice 
r G [2, 10]. Note a large proportion of parameter space has 
been excluded. 



Under the proper choice of KST parameters, shown 
shaded in Fig. one of the three types of ray (s = ±1, 0) 
may admit some prototypical coherent wave packet of po- 
larization \x, whose energy amplifies by 10^^ [i.e. Asfi > 
10^2]^ The clear region in Fig. indicates KST param- 
eters for which we have not yet found a transient which 
amphfies by 10'^^. 

D. Generalizations of our method which could 
generate stronger constraints on KST parameters 

With our extremely conservative approach - eliminat- 
ing those formulations with wave-packet solutions which 
amplify by 10'^^ in the future domain of dependence - 
we have already eliminated a broad region of parame- 
ter space. By relaxing some of our very restrictive as- 
sumptions, we expect we could discard still more KST 
parameters: 

1. Lower amplification cutoff: Currently, we require 
an enormous amplification before we eliminate a 
formulation; relaxing the requirement on amplifi- 
cation excludes more systems. 

2. Consider more transients Currently, we compute 
the amplification of only a handful of transients; a 
consideration of other transients (for example, in 
the neighborhood of circular photon orbits in PC) 
may allow us to exclude additional parameters. 

3. Consider a larger region Currently, we limit atten- 
tion only to the future domain of dependence of the 



initial data slice. Certain rays, however, remain 
within the computational domain for far longer. 
For example, in the PC case, rays near the horizon 
remain in the domain for arbitrarily long;^'' even 
the slowly-infalling rays last substantially longer 
than the domain of dependence. Therefore, by con- 
sidering the amplification of transients over a longer 
interval, we will discover significantly greater am- 
plification and thus exclude a significantly broader 
class of formulations of Einstein's equations. 



4. Combine with boundary conditions Finally, if we 
determine how geometric-optics solutions interact 
with boundary conditions, we can generalize our 
approach and address the late-time stability proper- 
ties of the evolution equations - or, in other words, 
address the stability properties of the full initial- 
plus-boundary value problem. 



One must take care to use the rays near the horizon in a sensible 
fashion. While analytically the rays remain within the computa- 
tional domain for arbitrarily long times, one cannot expect wave 
packet solutions to be resolved and present in a numerical solu- 
tion for arbitrarily long: the code has a finite smallest resolved 
scale. In practice, one must remember that whatever amplifi- 
cation one computes must be realistically attainable by some 
numerical simulation of fixed (though perhaps high) resolution 
in the coordinates of interest. 
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IX. CONCLUSIONS 

In this paper, we have demonstrated that certain tran- 
sients (prototypical coherent wave packets) can be used 
to veto a significant range of proposed formulations of 
Einstein's equations. We have described in considerable 
pedagogical detail precisely how to construct expressions 
for (or estimates of) the growth rate of prototypical co- 
herent wave packets [i.e. Eq. H20() ]. verify those esti- 
mates, and employ them to veto proposed formulations 
of Einstein's equations. These expressions employ no free 
parameters or knowledge of the solution, aside from a 
choice of plausible rays to examine. Moreover, despite 
the sometimes exhaustive details provided in Sections El 
and IVII the key tool - the growth rate of prototypical 
coherent wave packets [Eq. (|20|l ] - is easy to apply, with 
little conceptual, notational, or computational overhead 
(see, for example, the brief Sec. IVID H and its applica- 
tion in Sec. IVIII|I . Whether they are used conservatively, 
as in this paper, or generalized along the lines suggested 
in Sec. IVIII DI fi.e. using more rays and larger fragments 
of spacetime), we believe these techniques will provide a 
useful way to bound the number of proposed formulations 
before further tests are conducted (for example by the 
more ambitious Lindblom-Scheel energy-norm method) 
to decide whether a given formulation can produce effec- 
tive simulations. 

While our the specific examples of analyses in this pa- 
per have employed linearizations of the field equations 
themselves, we could just as well linearize a POSH system 
representing evolution equations for the constraint fields 
[see, for example, KST Eqs. (2.40-2.43)]. The evolution 
equations for the constraints have been emphasized by 
many other authors as an probe of unphysical behavior. 
Since the general arguments of Sections |TlllIIII and II VI do 
not depend on the precise FOSHLS used, we can perform 
a calculation following the same patterns as (for example) 
Sec. IVIIII to discover ill-behaved formulations.^^ 

In this paper, we have also discovered curious proper- 
ties of modes trapped on the horizon of a Schwarzchild 
hole in PG coordinates (Sec. IVIII|) . Analytically, we 
would expect that, if any growth rate for modes trapped 
on the horizon were positive, then these modes should 
grow without bound and be present in the evolution 
at late times. Numerically, however, we know that no 
resolved wave packets can appear at late times: such 
solutions would have to initiate arbitrarily close to the 
horizon, inconsistent with resolved, finite-resolution ini- 
tial data. Still, marginally-resolved solutions of similar 
character could potentially behave in an implementation- 



The author expects no new information can be obtained by such 
an analysis. Moreover, because the constraint equations, when 
written in first-order form, involve many more variables than 
the field equations themselves, such an analysis should prove 
substantially more challenging. 



dependent fashion, seeding outgoing modes which then 
propagate and amplify into the domain for all time. We 
shall explore this possibility in a future paper. 

APPENDIX A: USEFUL IDENTITIES USED IN 
THE TEXT 

1. An alternative approach to the group velocity 

The eigenequation which defines the natural polariza- 
tion spaces Bj and their associated eigenvalues ujj, 

for each Vj,a G Bj (cf. Sec. Ill Ap . may be alternatively 
expressed as 

A''{x)paPj{x,p) = UJj{x,p)Pj{x,p) . (Al) 

for Pj the projection operator to the jth eigenspace of 
A'^Pa- Differentiate this expression relative to pb, then 
apply Pj from the left, to find 

P^^^'P^ = - KP^ ■ (A2) 

Equivalently, if collection of basis vectors for 

the jth eigenspace which are orthonormal relative to the 
inner product generated by S, 

iv,,^,SA''v,,^) = S^pV; . (A3) 

2. The no-rotation condition 

In Section III C 21 we claim we can always find a basis 
for Bj at each point in the neighborhood of a given ray 
which satisfies the no-rotation condition [Eq. (|10() ]: 

{vj^[^,S{dt+A''da)v,,p])=0 (A4) 

for all a, In this section, we demonstrate explicitly 
how to construct a basis which satisfies both the no- 
rotation condition and which remains orthonormal. 

If the right-hand side term in Eq. ljA4|l is not zero in 
the basis Wj,q, we attempt to choose a new basis 

Vj^a = Raa^j^a 

such that Eq. 1)10(1 is satisfied by the new basis and more- 
over such that the new basis is orthonormal. In partic- 
ular, if the no- rotation condition is satisfied in the new 
basis, then (because PjA°-Pj = V^°) we know 

= ^i?aai?/3/3(l'j,[a,^(9t+A"aa)t;j-/3]) (A5) 
a 
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[In the second term above, the antisymmetrization is over 
only the barred indicies a and /?.] On the other hand, if 
the basis is orthonormal, then Yla^&a^pa = ^a/3j ™" 
plying 

a 

[In the above, the operator is symmetrized over the indi- 
cies a and /?.] Therefore, combining the two, we conclude 
that if the new basis is both orthonormal and satisfies the 
no-rotation, the matrix R must satisfy the ordinary dif- 
ferential equation 

{dt + V;da)Ro.o. = - {Vj,[o., S{dt + A''da)v,,f,]) Rs0 

subject to initial data Raa = Saa- Solutions for R and 
thus Vj^a exist in the neighborhood of a ray. 



3. Reorganizing inner products for the polarization 
equation 

In this section, we describe how to rearrange matters 
so the last term in Eq. (|7c|l - namely. 



Vj.a, S{dt+ A''da)Vj,f3 



(A6) 



- has a simpler form. If we choose a basis that satisfies 
the no-rotation condition [Eq. I|10|) ]. the antisymmetric 
part of this matrix is zero. Further, we may express the 
symmetric part of this expression by using the relations 

dt{Vj^a,SVj^p) = (Vj^a^SdtVj^p) (A7) 
+ {Vj^f3, SdtVj^a) 
+ {Vj,a,idtS)Vj^f3) 
daiVj^a^SA^Vj^ji) = {Vj^a,SA°'daVj^[3) 

+ {vj.p,SA''dtv,,^) (A8) 

+ {Vj,^,{daSA'')Vj,j,) 

[where we have observed that SA°' and A" are symmetric 
relative to the canonical inner product] and the expres- 
sions 



{Vj,a,SVj^l3) = Sap 
{Vj,a,SA''v,,p) = 6af3V; 



(A9) 
(AlO) 



[i.e. orthogonality and Eq. (jA3|) ]. These relations tell us 
that, if the no-rotation condition is satisfied, 

{Vj.a, S{dt+A''da)Vj,p)^ (All) 
iSa^pdaVix, k{x)) - i {V,,a, {dtS + daSA-) V,,p) 

Our notation for the first term on the right side (i.e. the 
divergence of the group velocity) is chosen to emphasize 
that the derivative da acts on all the dependence on a; - 
in particular, on any variation of ka = da4> with x. 



APPENDIX B: DEMONSTRATING THAT THE 

RAY-OPTICS APPROACH PROVIDES 
HIGH-QUALITY APPROXIMATE SOLUTIONS 
TO THE FOSHLS 

For any fixed initial data, the ray-optics solution ob- 
tained in Sec. will break down at some point along 
each ray. In this section, we estimate how long a solution 
obtained by solving Eq. ^ can be trusted. 

Specifically, in Sec. IB II we express the FOSHLS [Eq. 
©I using alternative variables better-suited to describ- 
ing the geometric-optics solution. Next, in Sec. IB 21 we 
survey the various orders of magnitude that arise in the 
problem. Using those orders of magnitude, in Sec. IB 31 
we estimate the error in Eq. (j^J that occurs when a 
geometric optics solution is substituted for u (e.g. we 
estimate how close the norm of the left side is to zero). 
Finally, knowing how much error we make when using a 
geometric-optics solution, in Sec. IB 4l we estimate how er- 
rors involved in a geometric-optics approximation grow; 
we therefore discover how long we can trust a purely 
geometric-optics-based evolution. 



1. Review: Writing equations in terms of (j) and di^a 

In Equations lO and © we describe how to parame- 
terize the A^-dimensional state vector u by using N func- 
tions di a and one additional function (p. If we insert this 
substitution into the original FOSHLS [Eq. lO], then dot 
the results against each of the orthonormal basis vectors 
vi^a, we obtain the equations 

= i{vi,a,Su)[dt + Vi''da]4> 



+ 1 Vi^a. S{dt+A''da~F)u^ 

= idiAdt + Vtda]4> (Bl) 

m p 

+ ^^d^Jvi^a, S{dt+ A^da- F)v.^ji 

m p ^ 

In the above, we have observed that A"" is symmetric 
relative to the inner product generated by S and that Vj^a 
is an eigenvector of A'^da4> with eigenvalue iVj = V-'dai'- 
We can further reorganize this equation by pulling out 
all terms that involve di explicitly, and also by using Eq. 
(HH) to simplify [vua,SA°-VLp) = VfSap- 



= idiAdt + V{'da]<l>+[dt + V{'da]di^. 



(B2) 



+ 2^ dup ( vi^a, S{dt + A- da - F)vi^p 

Vm,p)dad 



m^l P 
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2. Natural scales used in order-of-magnitude 
estimates 

To make order-of-magnitude arguments regarding the 
solution, we need to understand how the natural length 
scales of the problem enter into it. 

Rather than complicate the order-of-magnitude calcu- 
lation unnecessarily, we shall for simplicity proceed as 
if there existed only one characteristic speed. In other 
words, we shall freely convert between space and time 
units by using the norm of A"", for example, we can in- 
terpret Tj-nl^l as a natural length scale. Finally, for 
brevity, we shall assume space and time units are chosen 
so \A\ ~ 1. 

Even with the above simplification, many natural 
scales arise in the problem, including the magnitude of 
F; the natural length and time scales on which F and 
A vary; and the length scale on which the initial data 
varies. Again, for simplicity we shall summarize all these 
scales by only two numbers: 

• "Length" scale [L) We define the natural "length" 
scale to be the natural time scale that enters on 
the right side of Eq. (f7c|) . To be explicit, L is the 
smaller of \d\/\A\\dad\ and l/\F\. 

• "Variation" scale (Tvary) The remaining scales do 
not arise directly in the equation. They affect the 
propagation of the wave packet only because they 
determine the rate at which terms in the equation 
are modulated as the wave packet propagates in 
space and time. We therefore call the smallest of 
the remaining scales the variation scale (rvary); its 
value is the smallest of the length and time scales 
on which A and F vary. 



3. Degree to which ray-optics solution satisfies the 
FOSHLS 

Using the above rough estimates {L and Tvary) to char- 
acterize the magnitude of terms both used and neglected, 
we find that geometric-optics solutions [Eq. {T))] very 
nearly satisfy the full FOSHLS [Eq. Q, alternatively 
expressed as Eq. (|B2|) ]. To be explicit, when we insert 
a geometric-optics solution which propagates in the jth 
polarization [i.e. a solution to Eq. {T))] into Eq. I|lj2|l . we 
find the following: 

= idi^^{uoi-Uj)^[dt + V{'da]di^o. (B3) 

13 ^ ^ 

+ SA'^V^,p)dad 

m.JS 

+ S{dt+A''da-F)Vrn,pj 
7n=il ;3 ^ ^ 



[Here, we have used Eq. (|7a|l and the definition of ujj to 
simplify the first term.] 

We have two circumstances: 

• When I — j, the first three terms precisely cancel 
[see Eqs. (|7af) and iTTcl) ]. The only terms remaining 
are of order dm,f3 m ^ j. 

• On the other hand, when I j, the first term does 
not cancel. Rather, it is large, because da4> is large 
(i.e. we are in the short-wavelength limit), and the 
ll>i are proportional to da(t>- 

For brevity, assume the eigenvalues of A" arc of 
comparable magnitude, so to order of magnitude 
~ ~ w ~ ujj — ui . We may then express the 
equation when I ^ j as 

o = oM,^)-f o(d,-„/i) 

The second terms will force the first terms, gener- 
ally, to be nonzero. 

From the second case, we know that when I ^ j, \di\ ^ 
0{\dj\/Lijj). Combining this result with the first equa- 
tion, we conclude that when we use our trial solution, 
we are ignoring terms of order \dj\/L'^Lu when I — j and 
terms of order and \dj\/L when I ^ j. 

4. Length of time ray-optics solution can be trusted 

To estimate the integrated effects the neglected terms 
have on the diagonal and off-diagonal polarization am- 
plitudes dj^a and di^a, respectively, we integrate the pre- 
vious equations. 

When I j, we have a DE of form 

-^di^l3 + idi^l^Aij + Oi\dj\/L) ^ 

where we neglect smaller terms, where d/ds = dt + 
V°'da represents the derivative along a characteristic, and 
where Auj — uji — loj lo. Since we limit attention 
to Lo very large (i.e. cjTvary S> 1), we may ignore the 
weak effects of any time variation of L and treat it as 
constant. Since \dj\ varies along the characteristic much 
more slowly than does exp(— iws), we find that after an 
affine length T, |c?;| will be of order 

\di\ ^ T\dsdj\/LoL ^ \d,\T/L^uj . (B4) 

(Here, I assume L'^u is suitably averaged, as L will vary 
as the path evolves.) Similarly, when I — j, we ignore 
terms of order \dj\/L'^uj. We have a DE of form 

-^dj^a + (known, real) + 0{\dj\/ L'^ui) = 

Therefore, integrating along an affine length T of the ray, 
we expect errors in the dja's of relative magnitude 

Camp = T/L^LJ (B5) 
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when eamp is sail. In both cases, we see the neglected dual vector basis to select specific components, we find 
terms will be smaller than |c?j| by magnitude eamp- 



If we are simulating a given system, with fixed initial 
data, we can only trust a solution out to time T ~ L^tu. 
However, for any compact region of any characteristic 
(i.e. for any fixed T), we can always choose uj sufficiently 
large so the relative errors eamp is arbitrarily small. 



1 



= [dt + v;da + -daVf - o,^, D 



(Cl) 



APPENDIX C: WHEN DO PCWPS EXIST? 

Rather than evolve general wave packets using the full 
geometric-optics equations, for simplicity in this paper 
we often restrict attention to prototypcial coherent wave 
packets (PCWPs), which - if they exist - vastly simplify 
the problem of evolution (cf. Sec. IIII C(l . In this ap- 
pendix, we try to clarify the conditions under which pro- 
totypical coherent wave packet solutions exist as exact or 
approximate solutions to the geometric-optics equations. 

We can better understand under what conditions pro- 
totypical coherent wave packets exist if we rewrite the 
polarization equation [Eq. Ullfl ] using the basis /j''^ [Eq. 
(|14|l ]. When we do so, we find PCWP are exact so- 
lutions only for special circumstances. However, when 
some eigenvalue of Oj is large, we find that PCWPs arise 
naturally as limits of arbitrary coherent wave packets. 



Sufficient conditions for PCWP to be exact 
solution 



By definition, a prototypical coherent wave packet so- 
lution associated with the polarization direction fj^^ ex- 
ists only if there is a solution - exact or approximate 
- to Eq. (|(yH) with all Dj^, = except for v = ji [i.e. 
Dj^fi 7^ 0]. In other words, for a complete collection of 
PCWP solutions to exist, one for each /i, the third term 
in Eq. IjCip must be diagonal, or zero. 

[In fact, for the examples addressed in this paper (i.e. 
in Sections IVl and IVTI for propagation on the light cone), 
the third term is in fact exactly zero.] 



PCWP as limit of arbitrary rapidly-growing 
coherent wave packet 



1. Rewriting polarization equation in the basis of 
eigenvectors of Oj 

Basis vectors and their components: The basis vectors 
/j^'' are defined above. Since we express the polarization 
equation in component form relative to some no-rotation 
basis, we also need notation for the components /j^-* of 
these basis vectors relative to the no-rotation basis: 



(m) 



Dual vectors and their components: The basis vectors 
/j^'' are not necessarily orthogonal. To facilitate compu- 
tations, we define a dual basis /■^'' such that 



(ff\ Sf. 



The dual basis vectors can be expressed in terms of com- 
ponents, denoted relative to the no-rotation basis. 

Explicitly rewriting polarization equation: Substituting 
in the expansion 



J 3 



^ D f 



into the polarization equation [Eq. Hll|l ]. then using the 



If the largest eigenvalue Oji, of Oj is particularly large 
compared to the third term, then generic solutions to 
the polarization equation [Eq. will converge to a 

state with Dji, Dj^^ for [i ^ v (i.e. w — In 
other words, if the largest eigenvalue of Oj is large, then 
generic wave packets will converge to the PCWP with 
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APPENDIX D: BOUNDING THE ENERGY 
GROWTH RATE 



In Sections iHl and iHIl we introduced wave-packet solu- 
tions as solutions which are localized in the neighborhood 
of a given ray. However, while we obtaind expressions for 
the growth rate of certain specialized wave packet solu- 
tions [Sec. IIV| . in the main text of this paper we never 
provided a strict bound on the growth rate of a wave 
packet. 

In fact, we can bound the instantaneous growth rate of 
a coherent wave packet [Eq. (|19|l ] by a quantity indepen- 
dent of the precise polarization state w of that packet: 



1 dE (w, SQw) 

p-TT < max,„eB, -^ jtt" 

E at [w, Sw) 



(Dl) 



As Q is symmetric relative to S, it has a spectrum of 
real eigenvalues, each associated with eigenspaces that 
are orthogonal relative to S. It follows that if Kg is the 
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largest eigenvalue of Q, 



1 dE 
E~dt 



< 



(D2) 



This procedure follows precisely the same outline as the 
energy-norm upper bound discussed in LS Eq. (2.17) and 
(2.18). 

This upper bound on the growth rate for all polar- 
izations propagating along a given ray can be used as a 
line-by-line replacement for the maximum growth rate of 
PCWPs [Eq. (|20|l ] in practical calculations to determine 
the largest amplification possible by a wave packet prop- 
agating in the future domain of dependence [e.g. Sections 
rvTlland rVlTTj 



APPENDIX E: RAYS OPTICS AND KST 
2-PARAMETER FORMULATION 

While KST introduce a very large family of symmet- 
ric hyperbolic systems, they emphasize (and limit their 
calculations to) a simple 2-parameter subset. This two 
parameter system has both physical characteristic speeds 
and a simple principal part (i.e. simple form for A"'). As 
a result, the algebra required for its ray-optics limit (i.e. 
computations of ujj, etc) proves particularly simple. 



1. Generally 

The KST system has as variables the tensors gab, Pab, 
Mkab defined over 3-space, for a total of 6-|-6-|-18=30 
fields. 



a. Principal part and symmetrizer 

The principal part has the form [KST Eq. (2.59), along 
with the definition of do in KST Eq. (2.10)]: 

i^t-|3''^a)g^J ^ (Ela) 

i^t-|3''^a)P^J+Ng-''^aMb^J ~ (Elb) 
(dt - P''da)Mk^, + N^kP^J ^ (Elc) 

After linearizing about a background solution, this prin- 
cipal part and a choice of representation for the fields (i.e. 
u) gives us the explicit form for A°'. We may represent 
the result as 



(E2) 



for 1 the identity operator and a construction which 
depends only on the background metric g and the choice 
of field ordering used in going to a matrix representation 
(i.e. the order of the fields in u). 

This principal part is symmetric hyperbolic, using as 
symmetrizer (for example) LS Eq. (3.67): 



(u, Su) = g'"'g'"'dgabdg-ab + g'^'^ g"' dPabdP^, 



bb^ 



This symmetrizer (represented here in LS notation) 
amounts to nothing more than the naturally-constructed 
sum of squares of components of 5, P, and M. 



b. Eigenvalues and group velocity 

From the principal part, we can deduce the three pos- 
sible eigenvalues: 



Pb 



(E4) 



where s = 0, ±1. From this expression we obtain the 
group velocities 



<(a;,p) = 
where = g'^^pb / \lW^p7Ps 



-(3" + sNp" 



(E5) 



c. Eigenfields and projection operators 

KST tabulate the eigenfields of the principal-part op- 
erator [Eq. (|Eip ] in KST Eq. (2.61) and the surrounding 
text. These expressions yield the following orthonormal 
basis vectors for the three eigenspaces of A'^pa : 



"^o.g.ab 
'^o^x^ab 

'^o.y ,ab 

V±,ab 



gab 

[MqabXl - PuX^p'^Mqab] 
v/1 - {p^Xa^ 

[Mqaby'i - PurP^Mgab] 
y/1 -{p'^Va? 



V2 



[Pab ± p''Mqab\ 



(E6a) 
(E6b) 

(E6c) 

(E6d) 



+9 



aa ^hb ^kk 



dMuahdM, 



kab 



(E3) 



where we treat symbols for the fields as basis vectors in 
the space of fields, and where x and y are two 3-vector 
fields not parallel to p and which are orthonormal relative 
to the metric gab at each point. 



2. Special case: Flat spatial metric 

When the KST system is applied to a time- 
independent solution with a flat spatial metric, the al- 
gebra simplifies substantially. For example, the sym- 
metrizer S [Eq. IjESp ] is the identity operator on the set 
of fields. The inner product generated by S is therefore 
constant in space and time. 



a. Simplifying general polarization equation 

Since we fix the basis vector convention by Eq. HE6p . 
we must use the polarization equation in the form of Eq. 
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(f7c)l (i.e. we do not generically expect the no-rotation 
condition to hold) . We therefore must evaluate the term 

{Vj,p, S{dt + A'^da)Vj,0) = {Vjj,, S{dt - f3''da)v,,0) 

+ N{Vj,p,SAldaV,,p) 

Since we know how the basis vectors change as a function 
of the congruence direction k — p [Eqs. IjEBallEBdji ]. we 
can rewrite this expression in terms of our knowledge of 
the congruence and (3°". 

For example, for the fields propagating forward along 
the congruence at unit speed (j = s = ±1), we have 

{Vs,ab, SA'^^daVs^cd) = (l^dgpA SacSbd (E7) 



(E8) 



identity, so we just differentiate the result we obtain from 
Eq. 



Now, if we take the inner product over the fields, we 
arrive at the expression 



{Vj,a, idaSA'')Vj^p) = Sap 



\P\ 



(E9) 
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